function [UCoef,ExtraInfo]...
        =Bellman_U(PP,SS,Pir,dE,T,Profit,TAU,ir_dom,ir_ext,w_N,w_H,Lag_PortfolioInfo,PortfolioInfo,UCoefp,Flag_SS)
% Note:
%   This Bellman equation maps from EVFunCoefp to EVFunCoef.
% Profit=SS.Profit;ir_dom=SS.ir_dom;ir_ext=SS.ir_ext;Pir=SS.Pir;dE=SS.dE;
% T=SS.T; w_N=SS.w_N; w_H=SS.w_H; State=SS.FunApp.V.State;TAU=SS.TAU;
% PortfolioInfo=SS.PortfolioInfo; Lag_PortfolioInfo=SS.PortfolioInfo;
%% Sequential Solver
% Step 1: from UCoefp to VbarCoef
[Opt_Vbar,Vbar]     =   Fun_Vbar(PP,SS,T,Profit,TAU,ir_dom,ir_ext,w_N,w_H,PortfolioInfo,UCoefp,SS.FunApp.V.State,Flag_SS);
VbarCoef_Mat        =   SS.FunApp.V.BasMat_0\reshape(Opt_Vbar,SS.FunApp.V.MatSize_State);
VbarCoef            =   VbarCoef_Mat(:);
% Step 2: from VbarCoef to UCoef
[Opt_U,U]           =   Fun_U(PP,SS,Pir,dE,Lag_PortfolioInfo,VbarCoef,SS.FunApp.U.State,Flag_SS);
UCoef_Mat           =   SS.FunApp.U.BasMat_0\reshape(Opt_U,SS.FunApp.U.MatSize_State);
UCoef               =   UCoef_Mat(:);
%% Collocation Operation
AccType     =   'Howard';
if Flag_SS==1 
    switch AccType
        case 'Collocation'
            % Collocation
            % VBasMat*VbarCoef = FlowVec+CoefMat*UCoefp;
            FlowVec     =   Vbar.FlowVec;
            CoefMat     =   Vbar.CoefMat;
            % from V to U
            FlowVec     =   U.FlowVec+U.CoefMat*( SS.FunApp.V.BasMat_0_State\FlowVec );
            CoefMat     =   U.CoefMat*( SS.FunApp.V.BasMat_0_State\CoefMat );
            
            UCoef       =   ( SS.FunApp.U.BasMat_0_State-CoefMat )\FlowVec;
        case 'Howard'
            % Howard
            HowardItNum     =   40;
            for ii=1:HowardItNum
                Temp_V      =   Vbar.FlowVec+Vbar.CoefMat*UCoefp;
                Temp_VbarCoef=  SS.FunApp.V.BasMat_0_State\Temp_V;
                [Temp_U,U]  =   Fun_U(PP,SS,Pir,dE,Lag_PortfolioInfo,Temp_VbarCoef,SS.FunApp.U.State,Flag_SS);
                UCoef_Mat   =   SS.FunApp.U.BasMat_0\reshape(Temp_U,SS.FunApp.U.MatSize_State);
                UCoef       =   UCoef_Mat(:);

                Diff        =   max(abs(UCoef-UCoefp));
                if Diff<1e-8
                    break;
                end
                UCoefp     =   UCoef;
            end
    end
end
%% Collect the Information
ExtraInfo   =   struct('FunCoef',struct('Vbar',VbarCoef,'U',UCoef),...
                       'Vbar',Vbar,'U',U);
